import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import igrins_mod as ig # Custom module file for igrins shenanigans
import glob
import os
import warnings
import astropy.units as u
from astroquery.nist import Nist # atomic lines
# from astroquery.linelists.cdms import CDMS # molecular lines?
from lmfit import Model, Parameters, CompositeModel
from lmfit.models import GaussianModel
from scipy.stats import chisquare
from scipy.integrate import trapz, simpson
from scipy.optimize import curve_fit
from astropy.io import fits
# Plotting Parameters
plt.rcParams['font.size'] = 20
plt.rcParams['axes.labelsize'] = 18
plt.rcParams['xtick.labelsize'] = 18
plt.rcParams['ytick.labelsize'] =18
plt.rcParams['legend.fontsize'] = 16
plt.rcParams['figure.titlesize'] = 20
plt.rcParams['axes.labelweight']='bold'
plt.rcParams['axes.linewidth'] = 3
plt.rcParams['xtick.major.size'] = 10
plt.rcParams['xtick.minor.visible'] = True
plt.rcParams['xtick.minor.size'] = 5
plt.rcParams['ytick.major.size'] = 10
plt.rcParams['ytick.minor.visible'] = True
plt.rcParams['ytick.minor.size'] = 5
plt.rcParams['ytick.direction'] = 'in'
plt.rcParams['xtick.direction'] = 'in'
%matplotlib inline
# Size of 1 spectral resolution element
# IGRINS Spectral Resolution
spec_res = 0.00001
c = 299792458 # speed of light m/s
# Use Normalized (single) Gaussian Distribution
def gaussian_func(x, ampl, center, std):
return ((ampl) / (std * np.sqrt(2 * np.pi)) * np.exp(-0.5 * ((x - center) / std) ** 2)) + 1
# Reduced and order-merged data filepath
# Laptop Path
data_path = "C:\\Users\\Savio\\Documents\\GitHub\\IGRINS-Spectra\\IGRINS_Merged"
# File path for figures to live in
fig_path = "C:\\Users\\Savio\\Documents\\GitHub\\IGRINS-SpectraIGRINS_figs\\standards_spectra"
# Create the folder if it doesn't exist
if not os.path.exists(fig_path):
os.makedirs(fig_path)
# Nicole's merged K-band spectra of some Taurus Standards
merged_standard_files = glob.glob(data_path + "/merged_standards/m*.fits")
standard_table = pd.read_csv('./standard_table_v2.txt', index_col=0) # csv of standards with file and Spectral Type, c/v TBA
# early_k = ('1', '2', '3', '4', '5')
# Only want to look at the K types for the moment
# standard_table_k = standard_table[standard_table['Spectral_Type'].str.startswith(('K1', 'K2', 'K5'))]
standard_shift = standard_table['Wavelength Shift'].values
standard_list = standard_table['Source'].values
hops_table = pd.read_csv('./hops_table.txt')
hops_list = hops_table['Source']
# Determine the maximum length of flux arrays for the standards
max_flux_length = max(len(fits.getdata(file)[1]) for file in standard_list)
max_wavelen_length = max(len(fits.getdata(file)[0]) for file in standard_list)
max_snr_length = max(len(fits.getdata(file)[2]) for file in standard_list)
hops_flux_length = max(len(fits.getdata(file)[1]) for file in hops_table['Source'])
hops_wavelen_length = max(len(fits.getdata(file)[0]) for file in hops_table['Source'])
hops_snr_length = max(len(fits.getdata(file)[2]) for file in hops_table['Source'])
# Initialize stacks with NaN values
wavelen_stack = np.full((max_wavelen_length, len(standard_list)), np.nan)
raw_flux_stack = np.full((max_flux_length, len(standard_list)), np.nan)
snr_stack = np.full((max_snr_length, len(standard_list)), np.nan)
raw_flux_err_stack = np.full((max_snr_length, len(standard_list)), np.nan)
# hops_wavelen_stack = np.full((max_wavelen_length, len(hops_table)), np.nan)
# hops_raw_flux_stack = np.full((max_flux_length, len(hops_table)), np.nan)
# hops_snr_stack = np.full((max_snr_length, len(hops_table)), np.nan)
# hops_raw_flux_err_stack = np.full((max_snr_length, len(hops_table)), np.nan)
# Fill stacks with data
for i, file in enumerate(standard_list):
# Get wavelength, flux, snr per resolution element data
data = fits.getdata(file)
wavelen, flux, snr = data[0], data[1], data[2]
# Clean data a bit
snr_min = 50 # Minimum SNR
snr_max = 1e4 # Maximum SNR
snr_cut = (snr > snr_min) # & (snr < snr_max) # Bitwise SNR masking
flux_min = 10000 # Minimum flux
flux_cut = flux > flux_min # Bitwise flux masking
wavelen_min = 2.0
wavelen_max = 2.4
wavelen_cut = (wavelen > wavelen_min) & (wavelen < wavelen_max)
# Apply masks and remove NaNs and infs
mask = flux_cut & wavelen_cut & snr_cut
wavelen = wavelen[mask]
flux = flux[mask]
snr = snr[mask]
# Remove NaNs and infs from wavelen, flux, and snr arrays
valid_indices = ~np.isnan(wavelen) & ~np.isnan(flux) & ~np.isnan(snr)
wavelen = wavelen[valid_indices] - standard_shift[i]
flux = flux[valid_indices]
snr = snr[valid_indices]
# Check for NaNs in the final arrays
if np.any(np.isnan(wavelen)) or np.any(np.isnan(flux)) or np.any(np.isnan(snr)):
print(f"NaNs found in data for file {file} after cleaning")
wavelen_stack[:len(wavelen), i] = wavelen # Wavelength arrays for each standard
raw_flux_stack[:len(flux), i] = flux
snr_stack[:len(snr), i] = snr
raw_flux_err_stack[:len(flux), i] = flux / snr
# for i, file in enumerate(hops_list):
# # Get wavelength, flux, snr per resolution element data
# data = fits.getdata(file)
# wavelen, flux, snr = data[0], data[1], data[2]
# # Clean data a bit
# snr_min = 50 # Minimum SNR
# snr_max = 1e4 # Maximum SNR
# snr_cut = (snr > snr_min) # & (snr < snr_max) # Bitwise SNR masking
# flux_min = 10000 # Minimum flux
# flux_cut = flux > flux_min # Bitwise flux masking
# wavelen_min = 2.0
# wavelen_max = 2.4
# wavelen_cut = (wavelen > wavelen_min) & (wavelen < wavelen_max)
# # Apply masks and remove NaNs and infs
# mask = flux_cut & wavelen_cut & snr_cut
# wavelen = wavelen[mask]
# flux = flux[mask]
# snr = snr[mask]
# # Remove NaNs and infs from wavelen, flux, and snr arrays
# valid_indices = ~np.isnan(wavelen) & ~np.isnan(flux) & ~np.isnan(snr)
# wavelen = wavelen[valid_indices] - standard_shift[i]
# flux = flux[valid_indices]
# snr = snr[valid_indices]
# # Check for NaNs in the final arrays
# if np.any(np.isnan(wavelen)) or np.any(np.isnan(flux)) or np.any(np.isnan(snr)):
# print(f"NaNs found in data for file {file} after cleaning")
# hops_wavelen_stack[:len(wavelen), i] = wavelen # Wavelength arrays for each standard
# hops_raw_flux_stack[:len(flux), i] = flux
# hops_snr_stack[:len(snr), i] = snr
# hops_raw_flux_err_stack[:len(flux), i] = flux / snr
# Directly query NIST to find line features in K-band
with warnings.catch_warnings(): # Ignore warnings
warnings.simplefilter('ignore')
lines_table = Nist.query(2.08*u.um,2.35*u.um,
linename = 'Na I, Sc I, Si I, Fe I, Fe II, Al I, Mg I, Ca I, H I, Ti I',
energy_level_unit='eV',output_order='wavelength')
igrins_wav_cut = (lines_table['Observed'] > 2.08) & (lines_table['Observed'] < 2.35)
lines_table = lines_table[igrins_wav_cut]
# lines_table = pd.read_csv('lines_table.txt')
# Make masks for the table of all the lines just in case I want to peek at certain transitions/wavelengths
na1_mask = lines_table['Spectrum'] == 'Na I'
sc1_mask = lines_table['Spectrum'] == 'Sc I'
si1_mask = lines_table['Spectrum'] == 'Si I'
fe1_mask = lines_table['Spectrum'] == 'Fe I'
fe2_mask = lines_table['Spectrum'] == 'Fe II'
al1_mask = lines_table['Spectrum'] == 'Al I'
mg1_mask = lines_table['Spectrum'] == 'Mg I'
ca1_mask = lines_table['Spectrum'] == 'Ca I'
h1_mask = lines_table['Spectrum'] == 'H I'
ti1_mask = lines_table['Spectrum'] == 'Ti I'
# Just add all the masks to a list for the sake of my plotting a few cells down
mask_list = [na1_mask,sc1_mask,si1_mask,fe1_mask,al1_mask,mg1_mask,ca1_mask,h1_mask,ti1_mask]
color_list = ['purple', 'orange', 'green', 'blue', 'brown', 'crimson', 'olive', 'cyan', 'darkgreen']
Fe I ~ 2.084 (lines_table[fe1_mask][6]), 2.275 (lines_table[fe1_mask][109]), 2.284 (lines_table[fe1_mask][115])
Si I ~ 2.092 (lines_table[si1_mask][0])
Mg I ~ 2.106 (lines_table[mg1_mask][0]), 2.282 lines_table[mg1_mask][11]
Al I ~ 2.11 (lines_table[al1_mask][1])
Ca I ~ 2.283
Ti I ~ 2.29 (lines_table[ti1_mask][69])
CO (2-0) > 2.294
lines_table[fe1_mask][6]
| Spectrum | Observed | Ritz | Transition | Rel. | Aki | fik | Acc. | Ei Ek | Lower level | Upper level | Type | TP | Line |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| str5 | float64 | str18 | float64 | str17 | float64 | float64 | str3 | str40 | str45 | str59 | str5 | str17 | str20 |
| Fe I | 2.08465 | 2.0846493 | 4796.968 | 1200 | -- | -- | -- | 6.01718438 - 6.61193286 | 3d7.(4F).5s | e 3F | 3 | 3d7.(4F).5p | 3F* | 4 | -- | -- | L11631 |
line_name = lines_table[fe1_mask][0]['Spectrum']
line_center = lines_table[fe1_mask][6]['Observed']
# from igrins_mod import local_continuum_fit
continuum_stack = []
norm_flux_stack = []
reg_idx_stack = []
# regions [left point,width]
regions = [(-200,10), (200,10)]
# number of regions I use and subtract 1 since index starts at 0
n = len(regions)-1
for i in range(len(standard_table)):
fig = plt.figure(figsize=(15,3))
continuum, region_indices = ig.local_continuum_fit(wavelen_stack[:,i],
raw_flux_stack[:,i],
poly_order = 1,
line_center = line_center,
spec_res = spec_res,
regions = regions)
# append indices to list of indices
reg_idx_stack.append(region_indices)
# append to list of the local continuum arrays
continuum_stack.append(continuum)
# normalize flux by dividing the flux by the continuum and create a list
norm_flux = raw_flux_stack[region_indices[0][0]:region_indices[n][1],i]/continuum_stack[i]
norm_flux_stack.append(norm_flux)
plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], raw_flux_stack[region_indices[0][0]:region_indices[n][1],i])
plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], continuum_stack[i])
for j in range(len(region_indices)):
plt.axvspan(wavelen_stack[region_indices[n-j][0],i],wavelen_stack[region_indices[n-j][1],i],color='red',alpha=.2)
plt.title(rf"{line_name} {line_center}$\mu m$ {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")
plt.show()
amplitude1_stack = []
center1_stack = []
sigma_stack = []
pcov_stack = []
best_model_stack = []
flux_constant = np.linspace(0,-1,len(standard_table))
# Define initial parameters for Gaussian fitting
init_params = (0.1, lines_table[fe1_mask][6]['Observed'], spec_res)
# wavelen_stack[region_indices[0][0]:region_indices[2][1],i],raw_flux_stack[region_indices[0][0]:region_indices[2][1],i]
for i in range(len(standard_table)):
# Perform Gaussian fitting for the current source
popt, pcov, param_err, best_model = ig.model_fit(ig.gaussian,
wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
norm_flux_stack[i], raw_flux_err_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
init_params,
maxfev = 50000)
amplitude1_stack.append(popt[0])
center1_stack.append(popt[1])
sigma_stack.append(popt[2])
pcov_stack.append(pcov)
best_model_stack.append(best_model)
# for i in range(len(standard_table)):
# fwhm = 2*np.sqrt(2*np.log(2))*sigma1_stack[i]
# area = np.abs(amplitude_stack[i]*sigma_stack[i]*np.sqrt(2*np.pi))
# ew = np.trapz(1-best_model_stack[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
fe1_2_0846_fit = []
for i in range(len(standard_table)):
fe1_2_0846_fit.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],amplitude1_stack[i],center1_stack[i],sigma_stack[i]))
area_int = trapz(1-best_model_stack[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
area_v2 = np.abs(amplitude1_stack[i])*sigma_stack[i]*np.sqrt(2*np.pi)
# print(f'{area_int} {area_v2}')
ew1_stack = [] # empty list to load in equivalent widths
# Integate 1-best_model to get the areas (equivalent widths) of each component Gaussian
for i in range(len(standard_table)):
fig = plt.figure(figsize=(15,3))
plt.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],norm_flux_stack[i], c='black')
plt.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],best_model_stack[i],c='red',lw=2,label='Best Model')
# plt.fill_between(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],best_model_stack[i],0,color='red',alpha=0.2)
# Fe I
# plt.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],fe1_2_0846_fit[i])#, label = 'Fe I 2.2263 Fit')
# plt.ylim(0)
plt.ylabel('Flux + Constant')
plt.xlabel(r'Wavelength [$\mu$m]')
plt.title(f"{line_name} ({line_center}$\mu$m) {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")
plt.legend(loc='best')
plt.show()
# Define the region for fitting
line_name = lines_table[fe1_mask][0]['Spectrum'] # Species
line_center = lines_table[fe1_mask][109]['Observed'] # Wavelength
# from igrins_mod import local_continuum_fit
continuum_stack = []
norm_flux_stack = []
reg_idx_stack = []
# regions [left point,width]
regions = [(-200,10), (200,10)]
# number of regions I use and subtract 1 since index starts at 0
n = len(regions)-1
for i in range(len(standard_table)):
fig = plt.figure(figsize=(15,3))
continuum, region_indices = ig.local_continuum_fit(wavelen_stack[:,i],
raw_flux_stack[:,i],
poly_order = 1,
line_center = line_center,
spec_res = spec_res,
regions = regions)
# append indices to list of indices
reg_idx_stack.append(region_indices)
# append to list of the local continuum arrays
continuum_stack.append(continuum)
# normalize flux by dividing the flux by the continuum and create a list
norm_flux = raw_flux_stack[region_indices[0][0]:region_indices[n][1],i]/continuum_stack[i]
norm_flux_stack.append(norm_flux)
plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], raw_flux_stack[region_indices[0][0]:region_indices[n][1],i])
plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], continuum_stack[i])
for j in range(len(region_indices)):
plt.axvspan(wavelen_stack[region_indices[n-j][0],i],wavelen_stack[region_indices[n-j][1],i],color='red',alpha=.2)
plt.title(f"{standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")
plt.show()
amplitude1_stack = []
center1_stack = []
sigma_stack = []
pcov_stack = []
best_model_stack = []
flux_constant = np.linspace(0,-1,len(standard_list))
# Define initial parameters for Gaussian fitting
init_params = (0.1, lines_table[fe1_mask][6]['Observed'], spec_res)
# wavelen_stack[region_indices[0][0]:region_indices[2][1],i],raw_flux_stack[region_indices[0][0]:region_indices[2][1],i]
for i in range(len(standard_table)):
# Perform Gaussian fitting for the current source
popt, pcov, param_err, best_model = ig.model_fit(ig.gaussian,wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
norm_flux_stack[i], raw_flux_err_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
init_params,
maxfev=5000)
amplitude1_stack.append(popt[0])
center1_stack.append(popt[1])
sigma_stack.append(popt[2])
pcov_stack.append(pcov)
best_model_stack.append(best_model)
# for i in range(len(standard_table)):
# fwhm = 2*np.sqrt(2*np.log(2))*sigma1_stack[i]
# area = np.abs(amplitude_stack[i]*sigma_stack[i]*np.sqrt(2*np.pi))
# ew = np.trapz(1-best_model_stack[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
# fe1_2_0846_fit = []
# for i in range(len(standard_table)):
# fe1_2_0846_fit.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],amplitude1_stack[i],center1_stack[i],sigma_stack[i]))
# area_int = trapz(1-best_model_stack[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
# area_v2 = np.abs(amplitude1_stack[i])*sigma_stack[i]*np.sqrt(2*np.pi)
# # print(f'{area_int} {area_v2}')
# ew1_stack = [] # empty list to load in equivalent widths
# Integate 1-best_model to get the areas (equivalent widths) of each component Gaussian
for i in range(len(standard_table)):
fig = plt.figure(figsize=(15,3))
plt.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],norm_flux_stack[i], c='black')
plt.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],best_model_stack[i],c='red',lw=2,label='Best Model')
# plt.fill_between(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],best_model_stack[i],0,color='red',alpha=0.2)
# Fe I
# plt.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],fe1_2_0846_fit[i])#, label = 'Fe I 2.2263 Fit')
# plt.ylim(0)
plt.ylabel('Flux + Constant')
plt.xlabel(r'Wavelength [$\mu$m]')
plt.title(f"{line_name} ({line_center}$\mu$m) {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")
plt.legend(loc='best')
plt.show()
c:\Users\Savio\anaconda3\envs\muler_dev\lib\site-packages\scipy\optimize\_minpack_py.py:1010: OptimizeWarning: Covariance of the parameters could not be estimated
warnings.warn('Covariance of the parameters could not be estimated',
# Define the region for fitting
line_name = lines_table[si1_mask][0]['Spectrum'] # Species
line_center = lines_table[si1_mask][0]['Observed'] # Wavelength
# from igrins_mod import local_continuum_fit
continuum_stack = []
norm_flux_stack = []
reg_idx_stack = []
# regions [left point,width]
regions = [(-200,10), (200,10)]
# number of regions I use and subtract 1 since index starts at 0
n = len(regions)-1
for i in range(len(standard_table)):
fig = plt.figure(figsize=(15,3))
continuum, region_indices = ig.local_continuum_fit(wavelen_stack[:,i],
raw_flux_stack[:,i],
poly_order = 1,
line_center = line_center,
spec_res = spec_res,
regions = regions)
# append indices to list of indices
reg_idx_stack.append(region_indices)
# append to list of the local continuum arrays
continuum_stack.append(continuum)
# normalize flux by dividing the flux by the continuum and create a list
norm_flux = raw_flux_stack[region_indices[0][0]:region_indices[n][1],i]/continuum_stack[i]
norm_flux_stack.append(norm_flux)
plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], raw_flux_stack[region_indices[0][0]:region_indices[n][1],i])
plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], continuum_stack[i])
for j in range(len(region_indices)):
plt.axvspan(wavelen_stack[region_indices[n-j][0],i],wavelen_stack[region_indices[n-j][1],i],color='red',alpha=.2)
plt.title(f"{line_name} {line_center} {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")
plt.show()
# Define the region for fitting
line_name = lines_table[mg1_mask][0]['Spectrum'] # Species
line_center = lines_table[mg1_mask][0]['Observed'] # Wavelength
# from igrins_mod import local_continuum_fit
continuum_stack = []
norm_flux_stack = []
reg_idx_stack = []
# regions [left point,width]
regions = [(-200,10), (200,10)]
# number of regions I use and subtract 1 since index starts at 0
n = len(regions)-1
for i in range(len(standard_table)):
fig = plt.figure(figsize=(15,3))
continuum, region_indices = ig.local_continuum_fit(wavelen_stack[:,i],
raw_flux_stack[:,i],
poly_order = 1,
line_center = line_center,
spec_res = spec_res,
regions = regions)
# append indices to list of indices
reg_idx_stack.append(region_indices)
# append to list of the local continuum arrays
continuum_stack.append(continuum)
# normalize flux by dividing the flux by the continuum and create a list
norm_flux = raw_flux_stack[region_indices[0][0]:region_indices[n][1],i]/continuum_stack[i]
norm_flux_stack.append(norm_flux)
plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], raw_flux_stack[region_indices[0][0]:region_indices[n][1],i])
plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], continuum_stack[i])
for j in range(len(region_indices)):
plt.axvspan(wavelen_stack[region_indices[n-j][0],i],wavelen_stack[region_indices[n-j][1],i],color='red',alpha=.2)
plt.title(f"{standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")
plt.show()
# Define the region for fitting
line_name = lines_table[al1_mask][0]['Spectrum'] # Species
line_center = lines_table[al1_mask][1]['Observed'] # Wavelength
# from igrins_mod import local_continuum_fit
continuum_stack = []
norm_flux_stack = []
reg_idx_stack = []
# regions [left point,width]
regions = [(-200,10), (200,10)]
# number of regions I use and subtract 1 since index starts at 0
n = len(regions)-1
for i in range(len(standard_table)):
fig = plt.figure(figsize=(15,3))
continuum, region_indices = ig.local_continuum_fit(wavelen_stack[:,i],
raw_flux_stack[:,i],
poly_order = 1,
line_center = line_center,
spec_res = spec_res,
regions = regions)
# append indices to list of indices
reg_idx_stack.append(region_indices)
# append to list of the local continuum arrays
continuum_stack.append(continuum)
# normalize flux by dividing the flux by the continuum and create a list
norm_flux = raw_flux_stack[region_indices[0][0]:region_indices[n][1],i]/continuum_stack[i]
norm_flux_stack.append(norm_flux)
plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], raw_flux_stack[region_indices[0][0]:region_indices[n][1],i])
plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], continuum_stack[i])
for j in range(len(region_indices)):
plt.axvspan(wavelen_stack[region_indices[n-j][0],i],wavelen_stack[region_indices[n-j][1],i],color='red',alpha=.2)
plt.title(f"{standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")
plt.show()
# Define the region for fitting
line_name = lines_table[ti1_mask][0]['Spectrum'] # Species
line_center = lines_table[ti1_mask][69]['Observed'] # Wavelength
# from igrins_mod import local_continuum_fit
continuum_stack = []
norm_flux_stack = []
reg_idx_stack = []
# regions [left point,width]
regions = [(-200,10), (200,10)]
# number of regions I use and subtract 1 since index starts at 0
n = len(regions)-1
for i in range(len(standard_table)):
fig = plt.figure(figsize=(15,3))
continuum, region_indices = ig.local_continuum_fit(wavelen_stack[:,i],
raw_flux_stack[:,i],
poly_order = 1,
line_center = line_center,
spec_res = spec_res,
regions = regions)
# append indices to list of indices
reg_idx_stack.append(region_indices)
# append to list of the local continuum arrays
continuum_stack.append(continuum)
# normalize flux by dividing the flux by the continuum and create a list
norm_flux = raw_flux_stack[region_indices[0][0]:region_indices[n][1],i]/continuum_stack[i]
norm_flux_stack.append(norm_flux)
plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], raw_flux_stack[region_indices[0][0]:region_indices[n][1],i])
plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], continuum_stack[i])
for j in range(len(region_indices)):
plt.axvspan(wavelen_stack[region_indices[n-j][0],i],wavelen_stack[region_indices[n-j][1],i],color='red',alpha=.2)
plt.title(f"{standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")
plt.show()
Sc I (2.2058 $\mu m$ and 2.2071 $\mu m$): index 1 and 2
Si I (2.2068 $\mu m$)
Mg/Al (~2.11 $\mu m$)
Fe I (2.2205–2.2346 $\mu m$) Over Ti region: index 10 and 11
Ca I (2.2614, 2.2631, 2.2657 $\mu m$): Only three lines
Na I Doublet (2.2062, 2.2089) & Sc I (2.2058)
# Define the region for fitting
line_name = lines_table[na1_mask][0]['Spectrum'] # Species
line_center = lines_table[na1_mask][0]['Observed'] # Wavelength
# from igrins_mod import local_continuum_fit
continuum_stack = []
norm_flux_stack = []
reg_idx_stack = []
# regions [left point,width]
regions = [(-200,30),(-140,30),(150,40),(350,30)]
# number of regions I use and subtract 1 since index starts at 0
n = len(regions)-1
poly_deg = 3
for i in range(len(standard_table)):
fig = plt.figure(figsize=(15,3))
continuum, region_indices = ig.local_continuum_fit(wavelen_stack[:,i],
raw_flux_stack[:,i],
poly_order = poly_deg,
line_center = line_center,
spec_res = spec_res,
regions = regions)
# append indices to list of indices
reg_idx_stack.append(region_indices)
# append to list of the local continuum arrays
continuum_stack.append(continuum)
# normalize flux by dividing the flux by the continuum and create a list
norm_flux = raw_flux_stack[region_indices[0][0]:region_indices[n][1],i]/continuum_stack[i]
norm_flux_stack.append(norm_flux)
plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], raw_flux_stack[region_indices[0][0]:region_indices[n][1],i])
plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], continuum_stack[i])
for j in range(len(region_indices)):
plt.axvspan(wavelen_stack[region_indices[n-j][0],i],wavelen_stack[region_indices[n-j][1],i],color='red',alpha=.2)
plt.title(f"Na Interval {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")
plt.show()
# Initialize storage lists
params_stack = [] # list of best of fit parameters
params_error_stack = [] # list of errors on each parameter
best_model_stack = [] # list of best fits to the dats
result_stack = [] # lmfit results
dely_stack = [] # lmfit errors of the fits
for i in range(len(standard_table)):
# Define initial parameters for Gaussian fitting
params = Parameters()
params.add('amp1', value=-1e-5, min=-1, max=0)
params.add('c1', value=lines_table[sc1_mask][19]['Observed'], min=lines_table[sc1_mask][19]['Observed']-0.5e-4,max=lines_table[sc1_mask][19]['Observed']+0.5e-4)
params.add('std1', value=spec_res)
params.add('amp2', value=-1e-5, min=-1, max=0)
params.add('c2', value=line_center, min=line_center-0.5e-4,max=line_center+0.5e-4)
params.add('std2', value=spec_res)
params.add('amp3', value=-1e-5, min=-1, max=0)
params.add('c3', value=lines_table[si1_mask][2]['Observed'],min=lines_table[si1_mask][2]['Observed']-0.5e-4,max=lines_table[si1_mask][2]['Observed']+0.5e-4)
params.add('std3', value=spec_res)
params.add('amp4', value=-1e-5, min=-2e-5, max=0)
params.add('c4', value=lines_table[sc1_mask][20]['Observed'],min=lines_table[sc1_mask][20]['Observed']-0.5e-4,max=lines_table[sc1_mask][20]['Observed']+0.5e-4)
params.add('std4', value=spec_res)
params.add('amp5', value=-1e-5, min=-1, max=0)
params.add('c5', value=lines_table[na1_mask][1]['Observed'],min=lines_table[na1_mask][1]['Observed']-0.5e-4,max=lines_table[na1_mask][1]['Observed']+0.5e-4)
params.add('std5', value=spec_res)
# mathematical/logic constraints on parameters
# enforce same widths for the Scandium and Sodium lines
# params['std4'].expr = 'std1'
# params['std5'].expr = 'std2'
# Adding a parameter for the wavelength differences
# params.add('delta', value=1e-3)
# Constraining the center parameters to have the same difference
# params['c2'].expr = 'c1 + delta'
# params['c3'].expr = 'c1 + 2*delta'
# params['c4'].expr = 'c1 + 3*delta'
# params['c5'].expr = 'c1 + 4*delta'
# Define the model
model = Model(ig.five_gaussian, nan_policy='omit')
# Perform the fit
result = model.fit(norm_flux_stack[i], params,
x=wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i],
weights=1 / raw_flux_err_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i])
result_stack.append(result) # lmfit ModelResult objects
dely = result.eval_uncertainty(sigma=5) # lmfit
dely_stack.append(dely)
# Append the best_model, error and parameters stacks
# could probably make a loop to the tune of for value in result etc...
params_stack.extend([result.params['amp1'].value, result.params['c1'].value, result.params['std1'].value,
result.params['amp2'].value, result.params['c2'].value, result.params['std2'].value,
result.params['amp3'].value, result.params['c3'].value, result.params['std3'].value,
result.params['amp4'].value, result.params['c4'].value, result.params['std4'].value,
result.params['amp5'].value, result.params['c5'].value, result.params['std5'].value])
params_error_stack.extend([result.params['amp1'].stderr, result.params['c1'].stderr, result.params['std1'].stderr,
result.params['amp2'].stderr, result.params['c2'].stderr, result.params['std2'].stderr,
result.params['amp3'].stderr, result.params['c3'].stderr, result.params['std3'].stderr,
result.params['amp4'].stderr, result.params['c4'].stderr, result.params['std4'].stderr,
result.params['amp5'].stderr, result.params['c5'].stderr, result.params['std5'].stderr])
best_model_stack.append(result.best_fit)
# print(result.fit_report())
params_arr = np.array(params_stack)
amps = params_arr[0::3]
amps_err = np.array(params_error_stack)[0::3]
centers = params_arr[1::3]
center_err = np.array(params_error_stack)[1::3]
sigmas = params_arr[2::3]
sigmas_err = np.array(params_error_stack)[2::3]
# Generate the fits for each Gaussian component
na1_2_2062_fit = []
sc1_2_2058_fit = []
si1_2_2068_fit = []
sc1_2_2071_fit = []
na1_2_2089_fit = []
num_gauss = 5
for i in range(len(standard_table)):
sc1_2_2058_fit.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], amps[0::num_gauss][i], centers[0::num_gauss][i], sigmas[0::num_gauss][i]))
na1_2_2062_fit.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], amps[1::num_gauss][i], centers[1::num_gauss][i], sigmas[1::num_gauss][i]))
si1_2_2068_fit.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], amps[2::num_gauss][i], centers[2::num_gauss][i], sigmas[2::num_gauss][i]))
sc1_2_2071_fit.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], amps[3::num_gauss][i], centers[3::num_gauss][i], sigmas[3::num_gauss][i]))
na1_2_2089_fit.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], amps[4::num_gauss][i], centers[4::num_gauss][i], sigmas[4::num_gauss][i]))
for i in range(len(standard_table)):
# Create subplots with adjusted size ratios
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(15, 6), gridspec_kw={'height_ratios': [4, 1]}, sharex=True)
ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], norm_flux_stack[i], c='black', alpha=0.2, label='Data')
ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], best_model_stack[i], c='red', lw=1, label='Best Model')
# ax1.errorbar(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], best_model_stack[i], yerr=dely_stack[i],fmt='ro', ms=3, lw=1, label='Best Model')
ax1.fill_between(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i],
best_model_stack[i]-dely_stack[i], best_model_stack[i]+dely_stack[i],
alpha=0.2, color='red')
# plot the locations of fitted centers
ax1.scatter(y=[1.03] * len(centers[i * num_gauss:(i + 1) * num_gauss]), x=centers[i * num_gauss:(i + 1) * num_gauss], marker="|", s=100, color='black')
ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], sc1_2_2058_fit[i], ls='-.',
label=r'$\bf{Sc~I~2.2058}$' "\n" rf"A={amps[0::num_gauss][i]:.3e}" "\n" rf"$\mu$={centers[0::num_gauss][i]:.5f}" "\n" rf"$\sigma$= {sigmas[0::num_gauss][i]:.3e}")
ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], na1_2_2062_fit[i], ls='-.',
label=r'$\bf{Na~I~2.2062}$' "\n" rf"A={amps[1::num_gauss][i]:.3e}" "\n" rf"$\mu$={centers[1::num_gauss][i]:.5f}" "\n" rf"$\sigma$= {sigmas[1::num_gauss][i]:.3e}")
ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], si1_2_2068_fit[i], ls='-.',
label=r'$\bf{Si~I~2.2068}$' "\n" rf"A={amps[2::num_gauss][i]:.3e}" "\n" rf"$\mu$={centers[2::num_gauss][i]:.5f}" "\n" rf"$\sigma$= {sigmas[2::num_gauss][i]:.3e}")
ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], sc1_2_2071_fit[i], ls='-.',
label=r'$\bf{Sc~I~2.2071}$' "\n" rf"A={amps[3::num_gauss][i]:.3e}" "\n" rf"$\mu$={centers[3::num_gauss][i]:.5f}" "\n" rf"$\sigma$= {sigmas[3::num_gauss][i]:.3e}")
ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], na1_2_2089_fit[i], ls='-.',
label=r'$\bf{Na~I~2.2089}$' "\n" rf"A={amps[4::num_gauss][i]:.3e}" "\n" rf"$\mu$={centers[4::num_gauss][i]:.5f}" "\n" rf"$\sigma$= {sigmas[4::num_gauss][i]:.3e}")
# ax1.set_ylim(top=1.05)
ax1.set_ylabel('Flux + Constant')
# Calculate residuals
residuals = norm_flux_stack[i] - best_model_stack[i]
sqsum_res = np.sum(residuals**2)
sqsum = np.sum((norm_flux_stack[i] - np.mean(norm_flux_stack[i]))**2)
R2 = 1 - (sqsum_res / sqsum)
# Plot residuals
ax2.scatter(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], residuals, c='red', label='Residuals', s=3)
ax2.axhline(y=0, color='grey', linestyle='--')
ax2.set_xlabel(r'Wavelength [$\mu$m]')
ax2.set_ylabel('res')
ax2.set_title(rf'$R^2$={R2:.4f}', fontsize=12)
plt.suptitle(rf"Na Interval {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")
ax1.legend(bbox_to_anchor=(1, 1), fontsize=12)
plt.show()
# print(result_stack[i].fit_report())
# params_stack = []
# params_error_stack = []
# pcov_stack = []
# best_model_stack = []
# for i in range(len(standard_table)):
# # Define initial parameters for Gaussian fitting
# init_params = (-0.05, lines_table[sc1_mask][19]['Observed'], spec_res,
# -0.5, line_center, spec_res,
# -0.05, lines_table[si1_mask][2]['Observed'], spec_res,
# -0.5, lines_table[na1_mask][1]['Observed'], spec_res)
# # bounds for each parameter
# # gauss_bounds = ([-1,-np.inf,-np.inf, -1,-np.inf,-np.inf, -1,-np.inf,-np.inf, -1,-np.inf,-np.inf],
# # [0,np.inf,np.inf, 0,np.inf,np.inf, 0,np.inf,np.inf, 0,np.inf,np.inf])
# # Perform Gaussian fitting for the current source
# # add error keyword then propogate for equiv width
# # compare equiv widths between objects of the same spectral types
# popt, pcov, param_err, best_model = ig.model_fit(ig.four_gaussian,
# wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
# norm_flux_stack[i],
# raw_flux_err_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
# init_params,
# maxfev=5000)
# # Append the best_model, error and parameters stacks
# params_error_stack.extend(param_err)
# params_stack.extend(popt)
# params_arr = np.array(params_stack)
# # best fit amplitudes
# amps = params_arr[0::3]
# amps_err = params_error_stack[0::3]
# # best fit centers(means)
# centers = params_arr[1::3]
# center_err = params_error_stack[1::3]
# # best fit widths
# sigmas = params_arr[2::3]
# sigmas_err = params_error_stack[2::3]
# pcov_stack.append(pcov)
# best_model_stack.append(best_model)
# na1_2_2062_fit = []
# sc1_2_2058_fit = []
# si1_2_2068_fit = []
# na1_2_2089_fit = []
# for i in range(len(standard_table)):
# sc1_2_2058_fit.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],amps[0::4][i],centers[0::4][i],sigmas[0::4][i]))
# na1_2_2062_fit.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],amps[1::4][i],centers[1::4][i],sigmas[1::4][i]))
# si1_2_2068_fit.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],amps[2::4][i],centers[2::4][i],sigmas[2::4][i]))
# na1_2_2089_fit.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],amps[3::4][i],centers[3::4][i],sigmas[3::4][i]))
# for i in range(len(standard_table)):
# # Create subplots with adjusted size ratios
# fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(15, 6), gridspec_kw={'height_ratios': [4, 1]}, sharex=True)
# ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], norm_flux_stack[i], c='black',alpha=0.2, label='Data')
# ax1.scatter(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], best_model_stack[i], color='red', label='Best Model')
# ax1.scatter(y=[1.03] * len(centers[i * 4:(i + 1) * 4]), x=centers[i * 4:(i + 1) * 4], marker="|", s=50, color='black')
# # for j in range(len(amps[::3])):
# ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], sc1_2_2058_fit[i],ls='-.',
# label= r'Sc I 2.2058 Fit' "\n" rf"A={amps[0::4][i]:.3e}" "\n" rf"$\mu$={centers[0::4][i]:.5f}" "\n" rf"$\sigma$= {sigmas[0::4][i]:.3e}")
# ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], na1_2_2062_fit[i],ls='-.',
# label= r'Na I 2.2062 Fit' "\n" rf"A={amps[1::4][i]:.3e}" "\n" rf"$\mu$={centers[1::4][i]:.5f}" "\n" rf"$\sigma$= {sigmas[1::4][i]:.3e}")
# ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], si1_2_2068_fit[i],ls='-.',
# label = r'Si I 2.2068 Fit' "\n" rf"A={amps[2::4][i]:.3e}" "\n" rf"$\mu$={centers[2::4][i]:.5f}" "\n" rf"$\sigma$= {sigmas[2::4][i]:.3e}")
# ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], na1_2_2089_fit[i],ls='-.',
# label= r'Na I 2.2089 Fit' "\n" rf"A={amps[3::4][i]:.3e}" "\n" rf"$\mu$={centers[3::4][i]:.5f}" "\n" rf"$\sigma$= {sigmas[3::4][i]:.3e}")
# ax1.set_ylim(top=1.05)
# ax1.set_ylabel('Flux + Constant')
# # Calculate residuals and R^2
# residuals = (norm_flux_stack[i] - best_model_stack[i])
# sqsum_res = np.sum(residuals**2)
# sqsum = np.sum((norm_flux_stack[i]-np.mean(norm_flux_stack[i])**2))
# R2 = 1-(sqsum_res/sqsum)
# # Plot residuals
# ax2.scatter(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], residuals, c='red', label='Residuals', s=3)
# ax2.axhline(y=0, color='grey', linestyle='--')
# # ax2.set_ylim(-0.1,0.1)
# ax2.set_xlabel(r'Wavelength [$\mu$m]')
# ax2.set_ylabel('res')
# ax2.set_title(rf'$R^2$={R2:.3f}', fontsize=12)
# plt.suptitle(rf"{line_name} ({line_center}$\mu$m) {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")
# ax1.legend(bbox_to_anchor=(1,1), fontsize=12)
# plt.show()
ew1_stack = [] # empty list to load in equivalent widths
ew2_stack = []
ew3_stack = []
ew4_stack = []
ew5_stack = []
# Integate 1-best_model to get the areas (equivalent widths) of each component Gaussian
for i in range(len(standard_table)):
# ew1 = (ig.gaussian_area(amps[0::5][i],sigmas[0::5][i])) # integrate the gaussian
# ew1_stack.append(ew1)
# ew2 = (ig.gaussian_area(amps[1::5][i],sigmas[1::5][i]))
# ew2_stack.append(ew2)
# ew3 = (ig.gaussian_area(amps[2::5][i],sigmas[2::5][i]))
# ew3_stack.append(ew3)
# ew4 = (ig.gaussian_area(amps[3::5][i],sigmas[3::5][i]))
# ew4_stack.append(ew4)
# ew5 = (ig.gaussian_area(amps[4::5][i],sigmas[4::5][i]))
# ew5_stack.append(ew5)
ew1 = np.trapz(1-na1_2_2062_fit[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i]) # integrate the gaussian
ew1_stack.append(ew1)
ew2 = np.trapz(1-sc1_2_2058_fit[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
ew2_stack.append(ew2)
ew3 = np.trapz(1-si1_2_2068_fit[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
ew3_stack.append(ew3)
ew4 = np.trapz(1-sc1_2_2071_fit[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
ew4_stack.append(ew4)
ew5 = np.trapz(1-na1_2_2089_fit[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
ew5_stack.append(ew5)
standard_table.loc[:,'ew_na1_2_2062'] = ew1_stack
standard_table.loc[:,'ew_sc1_2_2058'] = ew2_stack
standard_table.loc[:,'ew_si1_2_2068'] = ew3_stack
standard_table.loc[:,'ew_sc1_2_2071'] = ew4_stack
standard_table.loc[:,'ew_na1_2_2089'] = ew5_stack
fig = plt.figure(figsize=(15,5))
plt.plot(standard_table['Name'], standard_table['ew_sc1_2_2058'],ls='-',marker='o', label=f"Sc I 2.2058 $\mu$m")
plt.plot(standard_table['Name'], standard_table['ew_na1_2_2062'],ls='-',marker='o', label=f"Na I 2.2062 $\mu$m")
plt.plot(standard_table['Name'], standard_table['ew_si1_2_2068'],ls='-',marker='o', label=f"Si I 2.2068 $\mu$m")
plt.plot(standard_table['Name'], standard_table['ew_sc1_2_2071'],ls='-',marker='o', label=f"Sc I 2.2071 $\mu$m")
plt.plot(standard_table['Name'], standard_table['ew_na1_2_2089'],ls='-',marker='o', label=f"Na I 2.2089 $\mu$m")
# plt.scatter(standard_table['Spectral_Type'], area_stack)
# plt.ylim(np.nanmedian(ew_stack)+-0.001,np.nanmedian(ew_stack)+0.001)
# plt.ylim(1e-6,1e-3)
plt.xticks(rotation=45)
plt.xlabel('Spectral Type')
plt.ylabel(r'Equivalent Width [$\mu$m]')
# plt.title(f"{line_name} {line_center} $\mu$m")
# plt.yscale('log')
plt.legend(bbox_to_anchor=(1,1))
plt.show()
# Define the region for fitting
line_name = lines_table[mg1_mask][0]['Spectrum'] # Species
line_center = lines_table[mg1_mask][0]['Observed'] # Wavelength
# from igrins_mod import local_continuum_fit
continuum_stack = []
norm_flux_stack = []
reg_idx_stack = []
# regions [left point,width]
regions = [(-160,30),(-80,20),(100,20),(200,20),(270,20),(430,20)]
# number of regions I use and subtract 1 since index starts at 0
n = len(regions)-1
poly_deg = 3
for i in range(len(standard_table)):
fig = plt.figure(figsize=(15,3))
continuum, region_indices = ig.local_continuum_fit(wavelen_stack[:,i],
raw_flux_stack[:,i],
poly_order = poly_deg,
line_center = line_center,
spec_res = spec_res,
regions = regions)
# append indices to list of indices
reg_idx_stack.append(region_indices)
# append to list of the local continuum arrays
continuum_stack.append(continuum)
# normalize flux by dividing the flux by the continuum and create a list
norm_flux = raw_flux_stack[region_indices[0][0]:region_indices[n][1],i]/continuum_stack[i]
norm_flux_stack.append(norm_flux)
plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], raw_flux_stack[region_indices[0][0]:region_indices[n][1],i])
plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], continuum_stack[i])
for j in range(len(region_indices)):
plt.axvspan(wavelen_stack[region_indices[n-j][0],i],wavelen_stack[region_indices[n-j][1],i],color='red',alpha=.2)
plt.title(f"{line_name} {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")
plt.show()
plt.show()
# Initialize storage lists
params_stack = []
params_error_stack = []
best_model_stack = []
result_stack = []
dely_stack = []
# Define initial parameters for Gaussian fitting
init_params = (-0.5, line_center, spec_res,
-0.5, lines_table[al1_mask]['Observed'][0], spec_res)
for i in range(len(standard_table)):
# Define initial parameters for Gaussian fitting
params = Parameters()
params.add('amp1', value = 1e-5, max = 0)
params.add('c1', value = line_center)
params.add('std1', value = spec_res)
params.add('amp2', value = 1e-5, max = 0)
params.add('c2', value = lines_table[al1_mask]['Observed'][0])
params.add('std2', value = spec_res)
# params['std1'].expr = 'std2' # constrain widths to be the same
# Define Model
model = Model(ig.two_gaussian, nan_policy='omit')
result = model.fit(norm_flux_stack[i], params,
x = wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i],
weights = 1/raw_flux_err_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i])
result_stack.append(result) # lmfit ModelResult objects
dely = result.eval_uncertainty(sigma=5) # lmfit
dely_stack.append(dely)
# Append the best_model, error and parameters stacks
# could probably make a loop to the tune of for value in result etc...
params_stack.extend([result.params['amp1'].value, result.params['c1'].value, result.params['std1'].value,
result.params['amp2'].value, result.params['c2'].value, result.params['std2'].value])
params_error_stack.extend([result.params['amp1'].stderr, result.params['c1'].stderr, result.params['std1'].stderr,
result.params['amp2'].stderr, result.params['c2'].stderr, result.params['std2'].stderr])
best_model_stack.append(result.best_fit)
params_arr = np.array(params_stack)
amps = params_arr[0::3]
amps_err = params_error_stack[0::3]
centers = params_arr[1::3]
center_err = params_error_stack[1::3]
sigmas = params_arr[2::3]
sigmas_err = params_error_stack[2::3]
# Gaussian model fits for each source
mg1_2_1065_fits = []
al1_2_1098_fits = []
num_gauss = 2
# Getting the individual Gaussians into their own lists
for i in range(len(standard_table)):
mg1_2_1065_fits.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],amps[0::num_gauss][i],centers[0::num_gauss][i],sigmas[0::num_gauss][i]))
al1_2_1098_fits.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],amps[1::num_gauss][i],centers[1::num_gauss][i],sigmas[1::num_gauss][i]))
ew1_stack = []
ew2_stack = []
# ew3_stack = []
# ew4_stack = []
for i in range(len(standard_table)):
ew1 = np.trapz(1-mg1_2_1065_fits[i], wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
ew2 = np.trapz(1-al1_2_1098_fits[i], wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
ew1_stack.append(ew1)
ew2_stack.append(ew2)
standard_table['ew_mg1_2_1065'] = ew1_stack
standard_table['ew_al1_2_1098'] = ew2_stack
# Plotting each component Gaussian and the best model fit
for i in range(len(standard_table)):
# Create subplots with adjusted size ratios
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(15, 6), gridspec_kw={'height_ratios': [4, 1]}, sharex=True)
ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], best_model_stack[i], c='red', label='Gaussian Fit')
ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], norm_flux_stack[i], c='black', alpha=0.5)
ax1.fill_between(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
best_model_stack[i]-dely_stack[i], best_model_stack[i]+dely_stack[i],
alpha=0.2, color='red')
ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], mg1_2_1065_fits[i], ls=':',
label=r'Mg I 2.1065 $\mu$m' '\n' rf'A={amps[0::num_gauss][i]:.3e}' '\n' rf'$\mu$={centers[0::num_gauss][i]:.5f}' '\n' rf'$\sigma$={sigmas[0::num_gauss][i]:.3e}')
ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], al1_2_1098_fits[i], ls=':',
label=r'Al I 2.1098 $\mu$m' '\n' rf'A={amps[1::num_gauss][i]:.3e}' '\n' rf'$\mu$={centers[1::num_gauss][i]:.5f}' '\n' rf'$\sigma$={sigmas[1::num_gauss][i]:.3e}')
ax1.set_ylabel('Flux + Constant')
ax1.legend(bbox_to_anchor=(1,1))
# Calculate residuals
residuals = (norm_flux_stack[i] - best_model_stack[i])
sqsum_res = np.sum(residuals**2)
sqsum = np.sum((norm_flux_stack[i]-np.mean(norm_flux_stack[i])**2))
R2 = 1-(sqsum_res/sqsum)
# Plot residuals
ax2.scatter(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], residuals, c='red', label='Residuals', s=3)
ax2.axhline(y=0, color='grey', linestyle='--')
# ax2.set_ylim(-0.1,0.1)
ax2.set_xlabel(r'Wavelength [$\mu$m]')
ax2.set_ylabel('res')
ax2.set_title(rf'$R^2$={R2:.3f}', fontsize=12)
plt.suptitle(rf"Mg/Al Interval {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")
plt.show()
fig = plt.figure(figsize=(15,5))
plt.plot(standard_table['Spectral_Type'], standard_table['ew_mg1_2_1065'],ls='--',marker='p', label=f"Mg I 2.1065 $\mu$m")
plt.plot(standard_table['Spectral_Type'], standard_table['ew_al1_2_1098'],ls='--',marker='p', label=f"Al I 2.1098 $\mu$m")
# plt.plot(standard_table['Spectral_Type'], standard_table['ew_fe1_2_2838'],ls='--',marker='p', label=f"Fe I 2.2838$\mu$m")
# plt.scatter(standard_table['Spectral_Type'], area_stack)
# plt.ylim(np.nanmedian(ew_stack)+-0.001,np.nanmedian(ew_stack)+0.001)
# plt.ylim(-0.0001,0.0002)
# plt.yscale('log')
plt.xlabel('Spectral Type')
plt.ylabel(r'Equivalent Width [$\mu$m]')
# plt.title(f"{line_name} {line_center} $\mu$m")
plt.legend(bbox_to_anchor=(1,1))
plt.show()
# Define the region for fitting
line_name = lines_table[fe1_mask][0]['Spectrum'] # Species
line_center = lines_table[fe1_mask][85]['Observed'] # Wavelength
# from igrins_mod import local_continuum_fit
continuum_stack = []
norm_flux_stack = []
reg_idx_stack = []
# regions [left point,width]
regions = [(-600,30),(-350,50),(-150,60),(130,10),(400,100)]
# number of regions I use and subtract 1 since index starts at 0
n = len(regions)-1
poly_deg = 3
for i in range(len(standard_table)):
fig = plt.figure(figsize=(15,3))
continuum, region_indices = ig.local_continuum_fit(wavelen_stack[:,i],
raw_flux_stack[:,i],
poly_order = poly_deg,
line_center = line_center,
spec_res = spec_res,
regions = regions)
# append indices to list of indices
reg_idx_stack.append(region_indices)
# append to list of the local continuum arrays
continuum_stack.append(continuum)
# normalize flux by dividing the flux by the continuum and create a list
norm_flux = raw_flux_stack[region_indices[0][0]:region_indices[n][1],i]/continuum_stack[i]
norm_flux_stack.append(norm_flux)
plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], raw_flux_stack[region_indices[0][0]:region_indices[n][1],i])
plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], continuum_stack[i])
for j in range(len(region_indices)):
plt.axvspan(wavelen_stack[region_indices[n-j][0],i],wavelen_stack[region_indices[n-j][1],i],color='red',alpha=.2)
plt.title(f"{line_name} {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")
plt.show()
params_stack = []
params_error_stack = []
pcov_stack = []
best_model_stack = []
# Define initial parameters for Gaussian fitting
init_params = (-0.2, lines_table[ti1_mask]['Observed'][44], spec_res,
-0.2, lines_table[ti1_mask]['Observed'][45], spec_res,
-0.2, line_center, spec_res,
-0.1, lines_table[fe1_mask]['Observed'][86], spec_res,
-0.1, lines_table[ti1_mask]['Observed'][47], spec_res)
for i in range(len(standard_table)):
# Perform Gaussian fitting for the current source
popt, pcov, param_err, best_model = ig.model_fit(ig.five_gaussian,
wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
norm_flux_stack[i],
raw_flux_err_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
init_params,
maxfev=5000)
# Append the best_model, error and parameters stacks
params_error_stack.extend(param_err)
params_stack.extend(popt)
params_arr = np.array(params_stack)
amps = params_arr[0::3]
amps_err = params_error_stack[0::3]
centers = params_arr[1::3]
center_err = params_error_stack[1::3]
sigmas = params_arr[2::3]
sigmas_err = params_error_stack[2::3]
pcov_stack.append(pcov)
best_model_stack.append(best_model)
ti1_2_2217_fit = []
ti1_2_2239_fit = []
fe1_2_2263_fit = [] # maybe fix wavelength offsets
fe1_2_2266_fit = []
ti1_2_2280_fit = []
num_gauss = 5
for i in range(len(standard_table)):
ti1_2_2217_fit.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
amps[0::num_gauss][i],centers[0::num_gauss][i],sigmas[0::num_gauss][i]))
ti1_2_2239_fit.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
amps[1::num_gauss][i],centers[1::num_gauss][i],sigmas[1::num_gauss][i]))
fe1_2_2263_fit.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
amps[2::num_gauss][i],centers[2::num_gauss][i],sigmas[2::num_gauss][i]))
fe1_2_2266_fit.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
amps[3::num_gauss][i],centers[3::num_gauss][i],sigmas[3::num_gauss][i]))
ti1_2_2280_fit.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
amps[4::num_gauss][i],centers[4::num_gauss][i],sigmas[4::num_gauss][i]))
# Integate 1-best_model to get the areas (equivalent widths) of each component Gaussian
for i in range(len(standard_table)):
# Create subplots with adjusted size ratios
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(15, 6), gridspec_kw={'height_ratios': [4, 1]}, sharex=True)
ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], norm_flux_stack[i], alpha=0.2, c='black', label='Data')
ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], best_model_stack[i], c='red', lw=2, label='Best Model')
ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],ti1_2_2217_fit[i], ls='--',
label = 'Ti I 2.2217 Fit' '\n' rf'A={amps[0::num_gauss][i]:.3e}' '\n' rf'$\mu$={centers[0::num_gauss][i]:.5f}' '\n' rf'$\sigma$={sigmas[0::num_gauss][i]:.3e}')
ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],ti1_2_2239_fit[i], ls='--',
label = 'Ti I 2.2239 Fit' '\n' rf'A={amps[1::num_gauss][i]:.3e}' '\n' rf'$\mu$={centers[1::num_gauss][i]:.5f}' '\n' rf'$\sigma$={sigmas[1::num_gauss][i]:.3e}')
ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],fe1_2_2263_fit[i], ls='--',
label = 'Fe I 2.2263 Fit' '\n' rf'A={amps[2::num_gauss][i]:.3e}' '\n' rf'$\mu$={centers[2::num_gauss][i]:.5f}' '\n' rf'$\sigma$={sigmas[2::num_gauss][i]:.3e}')
ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],fe1_2_2266_fit[i], ls='--',
label = 'Fe I 2.2266 Fit' '\n' rf'A={amps[3::num_gauss][i]:.3e}' '\n' rf'$\mu$={centers[3::num_gauss][i]:.5f}' '\n' rf'$\sigma$={sigmas[3::num_gauss][i]:.3e}')
ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],ti1_2_2280_fit[i], ls='--',
label = 'Ti I 2.2280 Fit' '\n' rf'A={amps[4::num_gauss][i]:.3e}' '\n' rf'$\mu$={centers[4::num_gauss][i]:.5f}' '\n' rf'$\sigma$={sigmas[4::num_gauss][i]:.3e}')
ax1.set_ylabel('Flux + Constant')
ax1.legend(bbox_to_anchor=(1,1), fontsize=12)
# Calculate residuals
residuals = (norm_flux_stack[i] - best_model_stack[i])
sqsum_res = np.sum(residuals**2)
sqsum = np.sum((norm_flux_stack[i]-np.mean(norm_flux_stack[i])**2))
R2 = 1-(sqsum_res/sqsum)
# Plot residuals
ax2.scatter(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], residuals, c='red', label='Residuals', s=3)
ax2.axhline(y=0, color='grey', linestyle='--')
# ax2.set_ylim(-0.1,0.1)
ax2.set_xlabel(r'Wavelength [$\mu$m]')
ax2.set_ylabel('res')
ax2.set_title(rf'$R^2$={R2:.3f}', fontsize=12)
plt.suptitle(rf"Ti Interval {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")
plt.show()
ew1_stack = [] # empty list to load in equivalent widths
ew2_stack = []
ew3_stack = []
# Integate 1-best_model to get the areas (equivalent widths) of each component Gaussian
for i in range(len(standard_table)):
ew1 = np.trapz(1-fe1_2_2263_fit[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i]) # integrate the gaussian
ew1_stack.append(ew1)
ew2 = np.trapz(1-fe1_2_2266_fit[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
ew2_stack.append(ew2)
ew3 = np.trapz(1-ti1_2_2280_fit[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
ew3_stack.append(ew3)
standard_table.loc[:,'ew_fe1_2_2263'] = ew1_stack
standard_table.loc[:,'ew_fe1_2_2266'] = ew2_stack
standard_table.loc[:,'ew_ti1_2_2280'] = ew3_stack
fig = plt.figure(figsize=(15,5))
plt.plot(standard_table['Spectral_Type'], np.abs(standard_table['ew_fe1_2_2263']),ls='-',marker='o', label=f"Fe I 2.2263 $\mu$m")
plt.plot(standard_table['Spectral_Type'], np.abs(standard_table['ew_fe1_2_2266']),ls='-',marker='o', label=f"Fe I 2.2266 $\mu$m")
plt.plot(standard_table['Spectral_Type'], np.abs(standard_table['ew_ti1_2_2280']),ls='-',marker='o', label=f"Fe I 2.2280 $\mu$m")
plt.xlabel('Spectral Type')
plt.ylabel(r'Equivalent Width [$\mu$m]')
# plt.title(f"{line_name} {line_center} $\mu$m")
plt.legend(bbox_to_anchor=(1, 1))
plt.show()
# import matplotlib
# matplotlib.use('TkAgg') # Set the backend to TkAgg for interactivity
# # Define the region for fitting
# line_name = lines_table[ca1_mask][0]['Spectrum'] # Species
# line_center = lines_table[ca1_mask][0]['Observed'] # Wavelength
# # from igrins_mod import local_continuum_fit
# continuum_stack = []
# norm_flux_stack = []
# reg_idx_stack = []
# # Helper function to select regions interactively
# def select_regions(wavelengths, flux, num_regions):
# plt.plot(wavelengths, flux)
# plt.title('Click to select points for region exclusion')
# plt.xlabel('Wavelength')
# plt.ylabel('Flux')
# # Select points interactively
# selected_points = plt.ginput(n=num_regions*2, timeout=0) # n=num_regions*2 because we need start and end points for each region
# plt.close()
# regions = []
# for i in range(0, len(selected_points), 2):
# start = selected_points[i][0]
# end = selected_points[i+1][0]
# regions.append((start, end))
# return regions
# # number of regions to select interactively
# num_regions = 4
# poly_deg = 3
# for i in range(len(standard_table)):
# idx1 = np.searchsorted(wavelen_stack[:, i], line_center ) - 150
# idx2 = np.searchsorted(wavelen_stack[:, i], lines_table[ca1_mask]['Observed'][2]) + 150
# wavelengths = wavelen_stack[idx1:idx2, i]
# flux = raw_flux_stack[idx1:idx2, i]
# # Select regions interactively
# regions = select_regions(wavelengths, flux, num_regions)
# # Prepare the regions for the local_continuum_fit function
# regions = [(int(np.searchsorted(wavelengths, start)), int(np.searchsorted(wavelengths, end))) for start, end in regions]
# fig = plt.figure(figsize=(15, 3))
# continuum, region_indices = ig.local_continuum_fit(wavelengths,
# flux,
# poly_order=poly_deg,
# line_center=line_center,
# spec_res=spec_res,
# regions=regions)
# # Append indices to list of indices
# reg_idx_stack.append(region_indices)
# # Append to list of the local continuum arrays
# continuum_stack.append(continuum)
# # Normalize flux by dividing the flux by the continuum and create a list
# norm_flux = flux[region_indices[0][0]:region_indices[-1][1]] / continuum_stack[i]
# norm_flux_stack.append(norm_flux)
# plt.plot(wavelengths[region_indices[0][0]:region_indices[-1][1]], flux[region_indices[0][0]:region_indices[-1][1]])
# plt.plot(wavelengths[region_indices[0][0]:region_indices[-1][1]], continuum_stack[i])
# for j in range(len(region_indices)):
# plt.axvspan(wavelengths[region_indices[j][0]], wavelengths[region_indices[j][1]], color='red', alpha=.2)
# plt.title(f"{line_name} {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")
# plt.show()
# Define the region for fitting
line_name = lines_table[ca1_mask][0]['Spectrum'] # Species
line_center = lines_table[ca1_mask][0]['Observed'] # Wavelength
# from igrins_mod import local_continuum_fit
continuum_stack = []
norm_flux_stack = []
reg_idx_stack = []
# regions [left point,width]
regions = [(-100,50), (50,30), (240,50),(550,10)]
# number of regions I use and subtract 1 since index starts at 0
n = len(regions)-1
poly_deg = 3
for i in range(len(standard_table)):
fig = plt.figure(figsize=(15,3))
continuum, region_indices = ig.local_continuum_fit(wavelen_stack[:,i],
raw_flux_stack[:,i],
poly_order = poly_deg,
line_center = line_center,
spec_res = spec_res,
regions = regions)
# append indices to list of indices
reg_idx_stack.append(region_indices)
# append to list of the local continuum arrays
continuum_stack.append(continuum)
# normalize flux by dividing the flux by the continuum and create a list
norm_flux = raw_flux_stack[region_indices[0][0]:region_indices[n][1],i]/continuum_stack[i]
norm_flux_stack.append(norm_flux)
plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], raw_flux_stack[region_indices[0][0]:region_indices[n][1],i])
plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], continuum_stack[i])
for j in range(len(region_indices)):
plt.axvspan(wavelen_stack[region_indices[n-j][0],i],wavelen_stack[region_indices[n-j][1],i],color='red',alpha=.2)
plt.title(f"{line_name} {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")
plt.show()
np.abs(lines_table[ca1_mask]['Observed'][0] - lines_table[ca1_mask]['Observed'][1])
# lines_table[ca1_mask]['Observed'][1]
# lines_table[ca1_mask]['Observed'][2]
# lines_table[fe1_mask]['Observed'][104]
0.0017000000000000348
# Initialize storage lists
params_stack = []
params_error_stack = []
best_model_stack = []
result_stack = []
dely_stack = []
# Define initial parameters for Gaussian fitting
# init_params = (-0.5, lines_table[ca1_mask]['Observed'][0], spec_res,
# -0.5, lines_table[ca1_mask]['Observed'][1], spec_res,
# -0.5, lines_table[ca1_mask]['Observed'][2], spec_res,
# -0.5, lines_table[fe1_mask]['Observed'][104], spec_res)
from scipy.signal import savgol_filter
for i in range(len(standard_table)):
# Define initial parameters for Gaussian fitting
params = Parameters()
params.add('amp1', value = 1e-5, max = 0)
params.add('c1', value = lines_table[ca1_mask]['Observed'][0])
params.add('std1', value = spec_res)
params.add('amp2', value = 1e-5, max = 0)
params.add('c2', value = lines_table[ca1_mask]['Observed'][1])
params.add('std2', value = spec_res)
params.add('amp3', value = 1e-5, max = 0)
params.add('c3', value = lines_table[ca1_mask]['Observed'][2])
params.add('std3', value = spec_res)
params.add('amp4', value = 1e-5, max = 0)
params.add('c4', value = lines_table[fe1_mask]['Observed'][104])
params.add('std4', value = spec_res)
params['std1'].expr = 'std2' # constrain widths to be the same
params['std1'].expr = 'std3'
# params['std1'].expr = 'std2'
# Define Model
model = Model(ig.four_gaussian, nan_policy='omit')
result = model.fit(savgol_filter(norm_flux_stack[i], window_length=11, polyorder=3), params,
x = wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i],
weights = 1/raw_flux_err_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i])
result_stack.append(result) # lmfit ModelResult objects
dely = result.eval_uncertainty(sigma=5) # lmfit
dely_stack.append(dely)
# Append the best_model, error and parameters stacks
# could probably make a loop to the tune of for value in result etc...
params_stack.extend([result.params['amp1'].value, result.params['c1'].value, result.params['std1'].value,
result.params['amp2'].value, result.params['c2'].value, result.params['std2'].value,
result.params['amp3'].value, result.params['c3'].value, result.params['std3'].value,
result.params['amp4'].value, result.params['c4'].value, result.params['std4'].value])
params_error_stack.extend([result.params['amp1'].stderr, result.params['c1'].stderr, result.params['std1'].stderr,
result.params['amp2'].stderr, result.params['c2'].stderr, result.params['std2'].stderr,
result.params['amp3'].stderr, result.params['c3'].stderr, result.params['std3'].stderr,
result.params['amp4'].stderr, result.params['c4'].stderr, result.params['std4'].stderr])
best_model_stack.append(result.best_fit)
params_arr = np.array(params_stack)
amps = params_arr[0::3]
amps_err = params_error_stack[0::3]
centers = params_arr[1::3]
center_err = params_error_stack[1::3]
sigmas = params_arr[2::3]
sigmas_err = params_error_stack[2::3]
# Gaussian model fits for each source
ca1_2_2614_fits = []
ca1_2_2631_fits = []
ca1_2_2657_fits = []
fe1_2_2626_fits = []
num_gauss = 4
# Getting the individual Gaussians into their own lists
for i in range(len(standard_table)):
ca1_2_2614_fits.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], amps[0::num_gauss][i], centers[0::num_gauss][i], sigmas[0::num_gauss][i]))
ca1_2_2631_fits.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], amps[1::num_gauss][i], centers[1::num_gauss][i], sigmas[1::num_gauss][i]))
ca1_2_2657_fits.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], amps[2::num_gauss][i], centers[2::num_gauss][i], sigmas[2::num_gauss][i]))
fe1_2_2626_fits.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], amps[3::num_gauss][i], centers[3::num_gauss][i], sigmas[3::num_gauss][i]))
# Plotting each component Gaussian and the best model fit
for i in range(len(standard_table)):
# Create subplots with adjusted size ratios
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(15, 6), gridspec_kw={'height_ratios': [4, 1]}, sharex=True)
ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], best_model_stack[i], c='red', label=r'Best Fit $\pm 5 \sigma$')
ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], savgol_filter(norm_flux_stack[i], window_length=11, polyorder=3), c='black', alpha=0.5)
ax1.fill_between(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
best_model_stack[i]-dely_stack[i], best_model_stack[i]+dely_stack[i],
alpha=0.2, color = 'red' )
ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], ca1_2_2614_fits[i], ls='-.',
label=r'Ca I 2.2614 $\mu$m' "\n" rf"A={amps[0::num_gauss][i]:.3e}" "\n" rf"$\mu$={centers[0::num_gauss][i]:.5f}" "\n" rf"$\sigma$= {sigmas[0::num_gauss][i]:.3e}")
ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], ca1_2_2631_fits[i], ls='-.',
label=r'Ca I 2.2631 $\mu$m' "\n" rf"A={amps[0::num_gauss][i]:.3e}" "\n" rf"$\mu$={centers[0::num_gauss][i]:.5f}" "\n" rf"$\sigma$= {sigmas[0::num_gauss][i]:.3e}")
ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], ca1_2_2657_fits[i], ls='-.',
label=r'Ca I 2.2657 $\mu$m' "\n" rf"A={amps[0::num_gauss][i]:.3e}" "\n" rf"$\mu$={centers[0::num_gauss][i]:.5f}" "\n" rf"$\sigma$= {sigmas[0::num_gauss][i]:.3e}")
ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], fe1_2_2626_fits[i], ls='-.', c='magenta',
label=r'Fe I 2.2626 $\mu$m' "\n" rf"A={amps[0::num_gauss][i]:.3e}" "\n" rf"$\mu$={centers[0::num_gauss][i]:.5f}" "\n" rf"$\sigma$= {sigmas[0::num_gauss][i]:.3e}")
# plt.axvline()
plt.suptitle(rf"Ca Region {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")
ax1.set_ylabel('Flux + Constant')
ax1.legend(bbox_to_anchor=(1,1))
# Calculate residuals
residuals = (norm_flux_stack[i] - best_model_stack[i])
sqsum_res = np.sum(residuals**2)
sqsum = np.sum((norm_flux_stack[i]-np.mean(norm_flux_stack[i])**2))
R2 = 1-(sqsum_res/sqsum)
# Plot residuals
ax2.scatter(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], residuals, c='red', label='Residuals', s=3)
ax2.axhline(y=0, color='grey', linestyle='--')
# ax2.set_ylim(-0.1,0.1)
ax2.set_xlabel(r'Wavelength [$\mu$m]')
ax2.set_ylabel('res')
ax2.set_title(rf'$R^2$={R2:.3f}', fontsize=12)
plt.show()
ew1_stack = []
ew2_stack = []
ew3_stack = []
ew4_stack = []
for i in range(len(standard_table)):
ew1 = np.trapz(1-ca1_2_2614_fits[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
ew2 = np.trapz(1-ca1_2_2631_fits[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
ew3 = np.trapz(1-ca1_2_2657_fits[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
ew4 = np.trapz(1-fe1_2_2626_fits[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
ew1_stack.append(ew1)
ew2_stack.append(ew2)
ew3_stack.append(ew3)
ew4_stack.append(ew4)
standard_table['ew_ca1_2_2614'] = ew1_stack
standard_table['ew_ca1_2_2631'] = ew2_stack
standard_table['ew_ca1_2_2657'] = ew3_stack
standard_table['ew_fe1_2_2626'] = ew4_stack
fig = plt.figure(figsize=(15,5))
# plt.plot(standard_table['Spectral_Type'], standard_table['ew_sc1_2_2058'],ls='--',marker='p', label=f"Sc I 2.2058 $\mu$m")
# plt.plot(standard_table['Spectral_Type'], standard_table['ew_na1_2_2062'],ls='--',marker='p', label=f"Na I 2.2062 $\mu$m")
# plt.plot(standard_table['Spectral_Type'], standard_table['ew_si1_2_2068'],ls='--',marker='p', label=f"Si I 2.2068 $\mu$m")
# plt.plot(standard_table['Spectral_Type'], standard_table['ew_na1_2_2089'],ls='--',marker='p', label=f"Na I 2.2089 $\mu$m")
# plt.plot(standard_table['Spectral_Type'], np.abs(standard_table['ew_fe1_2_2263']),ls='-.',marker='o', label=f"Fe I 2.2263 $\mu$m")
# plt.plot(standard_table['Spectral_Type'], standard_table['ew_fe1_2_2266'],ls='-.',marker='o', label=f"Fe I 2.2266 $\mu$m")
plt.plot(standard_table['Spectral_Type'], standard_table['ew_ca1_2_2614'],ls=':',marker='s', label=f"Ca I 2.2614 $\mu$m")
plt.plot(standard_table['Spectral_Type'], standard_table['ew_ca1_2_2631'],ls=':',marker='s', label=f"Ca I 2.2631 $\mu$m")
plt.plot(standard_table['Spectral_Type'], standard_table['ew_ca1_2_2657'],ls=':',marker='s', label=f"Ca I 2.2657 $\mu$m")
plt.plot(standard_table['Spectral_Type'], standard_table['ew_fe1_2_2626'],ls=':',marker='s', label=f"Fe I 2.2626 $\mu$m")
# plt.scatter(standard_table['Spectral_Type'], area_stack)
# plt.ylim(np.nanmedian(ew_stack)+-0.001,np.nanmedian(ew_stack)+0.001)
# plt.ylim(-0.0001,0.0005)
# plt.yscale('log')
plt.xlabel('Spectral Type')
plt.ylabel(r'Equivalent Width [$\mu$m]')
# plt.title(f"{line_name} {line_center} $\mu$m")
plt.legend(bbox_to_anchor=(1, 1))
plt.show()
Mg I 2.2813 = lines_table[mg1_mask][8]
Fe I 2.2838 = lines_table[fe1_mask][115] (possible 116)
Ca ? = 2.285?
# Define the region for fitting
line_name = lines_table[mg1_mask][0]['Spectrum'] # Species
line_center = lines_table[mg1_mask][8]['Observed'] # Wavelength
# from igrins_mod import local_continuum_fit
continuum_stack = []
norm_flux_stack = []
reg_idx_stack = []
# regions [left point,width]
regions = [(-150,10),(-100,10),(70,10),(300,10),(500,10)]
# number of regions I use and subtract 1 since index starts at 0
n = len(regions)-1
poly_deg = 3
for i in range(len(standard_table)):
fig = plt.figure(figsize=(15,3))
continuum, region_indices = ig.local_continuum_fit(wavelen_stack[:,i],
raw_flux_stack[:,i],
poly_order = poly_deg,
line_center = line_center,
spec_res = spec_res,
regions = regions)
# append indices to list of indices
reg_idx_stack.append(region_indices)
# append to list of the local continuum arrays
continuum_stack.append(continuum)
# normalize flux by dividing the flux by the continuum and create a list
norm_flux = raw_flux_stack[region_indices[0][0]:region_indices[n][1],i]/continuum_stack[i]
norm_flux_stack.append(norm_flux)
plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], raw_flux_stack[region_indices[0][0]:region_indices[n][1],i])
plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], continuum_stack[i])
for j in range(len(region_indices)):
plt.axvspan(wavelen_stack[region_indices[n-j][0],i],wavelen_stack[region_indices[n-j][1],i],color='red',alpha=.2)
plt.title(f"{line_name} {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")
plt.show()
# min_wave = line_center - (100 * spec_res)
# max_wave = line_center + (100 * spec_res)
# Optimal parameters that will be appended from scipy curve fit below
amplitude1_stack = []
center1_stack = []
sigma1_stack = []
amplitude2_stack = []
centeR2_stack = []
sigma2_stack = []
amplitude3_stack = []
center3_stack = []
sigma3_stack = []
pcov_stack = [] # Covariant matrices
best_model_stack = [] # Best model fits
flux_constant = np.linspace(0,-1,len(standard_list))
# Define initial parameters for Gaussian fitting
init_params = ( -0.5, line_center, spec_res,
-0.5, lines_table[fe1_mask]['Observed'][114], spec_res,
-0.5, 2.2825, spec_res)
for i in range(len(standard_table)):
# Perform Gaussian fitting for the current source
popt, pcov, param_err, best_model = ig.model_fit(ig.three_gaussian,
wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
norm_flux_stack[i], raw_flux_err_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
init_params,
maxfev=5000)
# Append the Gaussian parameters, covariant matrices, and model fits as we loop through
amplitude1_stack.append(popt[0])
center1_stack.append(popt[1])
sigma1_stack.append(popt[2])
amplitude2_stack.append(popt[3])
centeR2_stack.append(popt[4])
sigma2_stack.append(popt[5])
amplitude3_stack.append(popt[6])
center3_stack.append(popt[7])
sigma3_stack.append(popt[8])
# amplitude4_stack.append(popt[7])
# center4_stack.append(popt[8])
pcov_stack.append(pcov)
best_model_stack.append(best_model)
# Gaussian model fits for each source
mg1_2_2813_fits = []
fe1_2_2818_fits = []
ca1_2_2825_fits = []
# Getting the individual Gaussians into their own lists
for i in range(len(standard_table)):
mg1_2_2813_fits.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],amplitude1_stack[i],center1_stack[i],sigma1_stack[i]))
fe1_2_2818_fits.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],amplitude2_stack[i],centeR2_stack[i],sigma2_stack[i]))
ca1_2_2825_fits.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],amplitude3_stack[i],center3_stack[i],sigma3_stack[i]))
ew1_stack = []
ew2_stack = []
ew3_stack = []
# ew4_stack = []
for i in range(len(standard_table)):
ew1 = np.trapz(1-mg1_2_2813_fits[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
ew2 = np.trapz(1-fe1_2_2818_fits[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
ew3 = np.trapz(1-ca1_2_2825_fits[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
ew1_stack.append(ew1)
ew2_stack.append(ew2)
ew3_stack.append(ew3)
standard_table['ew_mg1_2_2813'] = ew1_stack
standard_table['ew_fe1_2_2818'] = ew2_stack
standard_table['ew_ca1_2_2825'] = ew3_stack
# Plotting each component Gaussian and the best model fit
for i in range(len(standard_table)):
# Create subplots with adjusted size ratios
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(15, 6), gridspec_kw={'height_ratios': [4, 1]}, sharex=True)
ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], best_model_stack[i], c='red',label='Gaussian Fit')
ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], norm_flux_stack[i], c='black', alpha=0.5)
ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], mg1_2_2813_fits[i],ls=':',label=r'Mg I 2.2813 $\mu$m')
ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], fe1_2_2818_fits[i],ls=':',label=r'Fe I 2.2818 $\mu$m')
ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], ca1_2_2825_fits[i],ls=':',label=r'Ca ? 2.2825 $\mu$m')
# plt.axvline()
ax1.suptitle(rf"Mg, Ca?, Fe {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")
ax1.set_ylabel('Flux + Constant')
ax1.set_xlabel(r'Wavelength [$\mu$m]')
ax1.legend(bbox_to_anchor=(1,1))
plt.show()
c:\Users\Savio\anaconda3\envs\muler_dev\lib\site-packages\scipy\optimize\_minpack_py.py:1010: OptimizeWarning: Covariance of the parameters could not be estimated
warnings.warn('Covariance of the parameters could not be estimated',
--------------------------------------------------------------------------- AttributeError Traceback (most recent call last) Cell In[40], line 95 92 ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], ca1_2_2825_fits[i],ls=':',label=r'Ca ? 2.2825 $\mu$m') 94 # plt.axvline() ---> 95 ax1.suptitle(rf"Mg, Ca?, Fe {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}") 97 ax1.set_ylabel('Flux + Constant') 98 ax1.set_xlabel(r'Wavelength [$\mu$m]') AttributeError: 'Axes' object has no attribute 'suptitle'
fig = plt.figure(figsize=(15,5))
plt.plot(standard_table['Spectral_Type'], standard_table['ew_mg1_2_2813'],ls='--',marker='p', label=f"Mg I 2.2813 $\mu$m")
plt.plot(standard_table['Spectral_Type'], standard_table['ew_fe1_2_2818'],ls='--',marker='p', label=f"Fe I 2.2818 $\mu$m")
plt.plot(standard_table['Spectral_Type'], standard_table['ew_ca1_2_2825'],ls='--',marker='p', label=f"Ca I 2.2825 $\mu$m")
# plt.scatter(standard_table['Spectral_Type'], area_stack)
# plt.ylim(np.nanmedian(ew_stack)+-0.001,np.nanmedian(ew_stack)+0.001)
plt.ylim(-0.0001,0.0002)
# plt.yscale('log')
plt.xlabel('Spectral Type')
plt.ylabel(r'Equivalent Width [$\mu$m]')
# plt.title(f"{line_name} {line_center} $\mu$m")
plt.legend(bbox_to_anchor=(1,1))
plt.show()
# standard_table.to_csv('standard_table_v2.txt')